////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#ifndef BASEUTILS_H
#define BASEUTILS_H

#include "headers.hh"

class BaseUtils
{
  //===================================================================
  //                                Functions
  //===================================================================
 public:
  typedef Double_t (*FFunc)(const Double_t *xx, const Double_t *par);

  static Double_t BinomialD(const Double_t xx, const Double_t nn, const Double_t pp);
  static Double_t BinomialD(const Double_t *xx, const Double_t *par);
  static Double_t ProtectedExp(const Double_t xx);
  static Double_t ProtectedPow(const Double_t xx, const Double_t yy);

  static Double_t ResolutionGaus(const Double_t *xx, const Double_t *par);

  static Double_t GausN(const Double_t *xx, const Double_t *par, const Int_t ntype, FFunc basefunc, const Int_t npartype);

  //===================================================================
  //                               Algorithm
  //===================================================================
 public:
  static Bool_t SolveEq(Double_t &x1, Double_t &x2, const Double_t K1[], const Double_t K2[]);

  static Double_t GetCorrelatedError(const Double_t v0, const Double_t v1, const Double_t a0, const Double_t a1, const Double_t a2);

  static void Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint=0);

  //===================================================================
  //                                Fitting
  //===================================================================
 public:
  static TMatrixDSym * GetCov(TMinuit * tmn, const Bool_t ifixed[], const Double_t errs[]);

  //
  static Int_t BinnedLikelihoodFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[]=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  //
  static Int_t BinnedFitHist(const Bool_t kml, const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);

  //
  static Int_t ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, Double_t *chi=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static void SetFitPrintLevel(const Int_t lev){ fgFitPrintLevel = lev; }
  static void SetTol(const Double_t tol){fTol = tol;}
  static void SetMaxIter(const Int_t iter){fMaxIter = iter;}

 private:
  //
  static Int_t BinnedFitLevel0(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix,Double_t *cov, const Double_t *lows, const Double_t *highs);

  //
  static Int_t FitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);

  //chi2
  static Double_t GetFCNChisquare(const Double_t pars[]);
  //
  static void MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);

  //binned likelihood
  static Double_t GetFCNBinnedLikelihood(const Double_t pars[]);
  //
  static void MinFCNBinnedLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);

  //=================

  static Int_t fgFitPrintLevel;  //print level

  static Int_t fgFitNStep; //n step
  static Int_t fgFitNData; //n data
  static Double_t *fgFitX; //data x
  static Double_t *fgFitY; //data y
  static Double_t *fgFitEx; //err x
  static Double_t *fgFitEy; //err y
  static FFunc fgFitFunc; //fit func
  static Double_t fTol;
  static Int_t fMaxIter;
  //===================================================================
  //                               Histograms
  //===================================================================
 public:
  static Double_t FindYield(const TH2D *hh, const Double_t target);
  static Double_t LowMomAnchor(const Double_t xx, const Double_t anc, const Int_t opt);
  static TH1D* Regenerate(TH2D *hh,  Double_t anchorN, const Int_t ancOpt, const Double_t thres);
  static TH1D * GetHyield(const TH2D *h0, const TString hn);
  static void SelectXbins(TH2D *hh, const Int_t bins[]);
  static TH2D * CutHist(TH2D *hh, const Double_t fx0, const Double_t fx1, const Double_t fy0, const Double_t fy1, const Bool_t krest=kFALSE);
  static Int_t RebinX(TH2D * hh, const Int_t xrebin, const Double_t xmid, Int_t outbins[]);

  static TH1D * GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax, const Bool_t kbinlogx=0);
  static TH2D* NormalHist(const TH2D *h0, const Double_t thres=0, const Bool_t kmax=kFALSE);
  static Double_t * GetAxisArray(TAxis * aa);

  static TGraphAsymmErrors * TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0=(Char_t*)"", Char_t * sex1=(Char_t*)"", Char_t * sey0=(Char_t*)"", Char_t * sey1=(Char_t*)"", const Bool_t kXLogToLinear=kFALSE, const Bool_t kprint=kFALSE);

 private:
  static TH1D* ToPDF(const TH1 *h0, const TString hn);
  static void BinLogX(TAxis *axis);

  //===================================================================
  //                           end
  //===================================================================
};

#endif
